﻿#pragma once
#include <vector>
#include <algorithm>
#include <iostream>
#include <cmath>

class Matrix;

class Matrixf
{
public:
	// 对角矩阵
	using DiagonalMatrix = std::vector<double>;
	// 判断是否为方阵
	static bool IsSquareMatrix(const Matrix &matrix);
	// 创建一个矩阵
	static Matrix CreateMatrix(int row, int col);
	// 上下水平翻转
	static void FhudMatrix(Matrix &matrix);
	// 左右水平翻转
	static void FvlrMatrix(Matrix &matrix);
	// 方阵转置
	static Matrix TransposeMatrix(const Matrix &matrix);
	// 矩阵相加
	static Matrix AddMatrix(const Matrix &matrixA,
							const Matrix &matrixB);
	// 矩阵相减
	static Matrix SubMatrix(const Matrix &matrixA,
							const Matrix &matrixB);
	// 矩阵相乘
	static Matrix MulMatrix(const Matrix &matrixA,
							const Matrix &matrixB);
	// 数与矩阵相乘
	static Matrix MulMatrix(const double num, const Matrix &matrix);
	// 矩阵求逆
	static Matrix InverseMatrix(const Matrix &matrix);
	// 通过对角矩阵的一维数组转换为矩阵二维数组
	static Matrix CreateDigonal(const DiagonalMatrix &diagonalMatrix);
};

class Matrix
{
	friend class Matrixf;

public:
	using MatrixContainerType = std::vector<std::vector<double>>;

	Matrix();

	Matrix(int row, int col);

	Matrix(const MatrixContainerType &matrixContainer);

	Matrix(const Matrix &matrix);

	Matrix(const Matrixf::DiagonalMatrix &diagonalMatrix);

	Matrix operator*(const Matrix &matrix) const;

	Matrix operator+(const Matrix &matrix) const;

	Matrix operator-(const Matrix &matrix) const;

	std::vector<double> &operator[](int i);

	const std::vector<double> &operator[](int i) const;

	const MatrixContainerType &GetContianer() const;

	// 矩阵的行
	int Row() const;

	// 矩阵的列
	int Col() const;

	// 该矩阵的逆
	Matrix &Inverse();

	// 该矩阵的转置
	Matrix &Transform();

	void Show();

	void Normalize();

private:
	MatrixContainerType _matrixContainer;
};

double epsilon = 0;

double roundToFourDecimals(double value)
{
	return std::floor(value * 10000 + 0.5) / 10000;
}

bool Matrixf::IsSquareMatrix(const Matrix &matrix)
{
	if (matrix.Row() != matrix.Col())
		return false;

	return true;
}

Matrix Matrixf::CreateMatrix(int row, int col)
{
	Matrix ans(row, col);
	return ans;
}

void Matrixf::FhudMatrix(Matrix &matrix)
{
	std::reverse(matrix._matrixContainer.begin(), matrix._matrixContainer.end());
}

void Matrixf::FvlrMatrix(Matrix &matrix)
{
	for (int i = 0; i < (int)matrix.Row(); i++)
	{
		std::reverse(matrix[i].begin(), matrix[i].end());
	}
}

Matrix Matrixf::TransposeMatrix(const Matrix &matrix)
{
	Matrix ans = CreateMatrix(matrix.Col(), matrix.Row());
	int h = ans.Row();
	int l = ans.Col();
	for (int i = 0; i < h; i++)
	{
		for (int j = 0; j < l; j++)
		{
			ans[i][j] = matrix[j][i];
		}
	}
	return ans;
}

Matrix Matrixf::AddMatrix(const Matrix &matrixA,
						  const Matrix &matrixB)
{
	int h = matrixA.Row();
	int l = matrixA.Col();
	Matrix ans;
	ans = CreateMatrix(h, l);

	for (int i = 0; i < h; i++)
	{
		for (int j = 0; j < l; j++)
		{
			ans[i][j] = matrixA[i][j] + matrixB[i][j];
			if (abs(ans[i][j]) < epsilon)
			{
				ans[i][j] = 0;
			}
		}
	}
	return ans;
}

Matrix Matrixf::SubMatrix(const Matrix &matrixA,
						  const Matrix &matrixB)
{
	int h = matrixA.Row();
	int l = matrixA.Col();
	Matrix ans;
	ans = CreateMatrix(h, l);

	for (int i = 0; i < h; i++)
	{
		for (int j = 0; j < l; j++)
		{
			ans[i][j] = matrixA[i][j] - matrixB[i][j];
			if (abs(ans[i][j]) < epsilon)
			{
				ans[i][j] = 0;
			}
		}
	}
	return ans;
}

Matrix Matrixf::MulMatrix(const Matrix &matrixA,
						  const Matrix &matrixB)
{
	int matrixAH = matrixA.Row();
	int matrixAL = matrixA.Col();
	int matrixBH = matrixB.Row();
	int matrixBL = matrixB.Col();

	Matrix ans = CreateMatrix(matrixAH, matrixBL);
	for (int i = 0; i < matrixAH; i++)
	{
		for (int j = 0; j < matrixBL; j++)
		{
			ans[i][j] = 0;
			for (int k = 0; k < matrixAL; k++)
			{
				ans[i][j] += matrixA[i][k] * matrixB[k][j];
			}
			if (abs(ans[i][j]) < epsilon)
			{
				ans[i][j] = 0.0;
			}
		}
	}
	return ans;
}

Matrix Matrixf::MulMatrix(const double num, const Matrix &matrix)
{
	int matrixH = matrix.Row();
	int matrixL = matrix.Col();
	Matrix ans = CreateMatrix(matrixH, matrixL);
	for (int i = 0; i < matrixH; i++)
	{
		for (int j = 0; j < matrixL; j++)
		{
			ans[i][j] = num * matrix[i][j];
		}
	}
	return ans;
}

Matrix Matrixf::InverseMatrix(const Matrix &matrix)
{
	Matrix A(matrix);

	int n = A.Row();
	auto inv_A = CreateMatrix(n, n);
	auto L = CreateMatrix(n, n);
	auto U = CreateMatrix(n, n);
	auto inv_L = CreateMatrix(n, n);
	auto inv_U = CreateMatrix(n, n);
	// LU分解
	// L矩阵对角元素为1
	for (int i = 0; i < n; i++)
	{
		L[i][i] = 1;
	}
	// U矩阵第一行
	for (int i = 0; i < n; i++)
	{
		U[0][i] = A[0][i];
	}
	// L矩阵第一列
	for (int i = 1; i < n; i++)
	{
		L[i][0] = 1.0 * A[i][0] / A[0][0];
	}

	// 计算LU上下三角
	for (int i = 1; i < n; i++)
	{
		// 计算U（i行j列）
		for (int j = i; j < n; j++)
		{
			double tem = 0;
			for (int k = 0; k < i; k++)
			{
				tem += L[i][k] * U[k][j];
			}
			U[i][j] = A[i][j] - tem;
			if (abs(U[i][j]) < epsilon)
			{
				U[i][j] = 0.0;
			}
		}
		// 计算L（j行i列）
		for (int j = i; j < n; j++)
		{
			double tem = 0;
			for (int k = 0; k < i; k++)
			{
				tem += L[j][k] * U[k][i];
			}
			L[j][i] = 1.0 * (A[j][i] - tem) / U[i][i];
			if (abs(L[i][j]) < epsilon)
			{
				L[i][j] = 0.0;
			}
		}
	}
	// L U剩余位置设为0
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			if (i > j)
			{
				U[i][j] = 0.0;
			}
			else if (i < j)
			{
				L[i][j] = 0.0;
			}
		}
	}
	// LU求逆
	// 求矩阵U的逆
	for (int i = 0; i < n; i++)
	{
		inv_U[i][i] = 1 / U[i][i]; // U对角元素的值，直接取倒数
		for (int k = i - 1; k >= 0; k--)
		{
			double s = 0;
			for (int j = k + 1; j <= i; j++)
			{
				s = s + U[k][j] * inv_U[j][i];
			}
			inv_U[k][i] = -s / U[k][k]; // 迭代计算，按列倒序依次得到每一个值，
			if (abs(inv_U[k][i]) < epsilon)
			{
				inv_U[k][i] = 0.0;
			}
		}
	}
	// 求矩阵L的逆
	for (int i = 0; i < n; i++)
	{
		inv_L[i][i] = 1; // L对角元素的值，直接取倒数，这里为1
		for (int k = i + 1; k < n; k++)
		{
			for (int j = i; j <= k - 1; j++)
			{
				inv_L[k][i] = inv_L[k][i] - L[k][j] * inv_L[j][i];
				if (abs(inv_L[k][i]) < epsilon)
				{
					inv_L[k][i] = 0.0;
				}
			}
		}
	}

	inv_A = MulMatrix(inv_U, inv_L);

	return inv_A;
}

Matrix Matrixf::CreateDigonal(const DiagonalMatrix &diagonalMatrix)
{
	Matrix ans = CreateMatrix(diagonalMatrix.size(), diagonalMatrix.size());

	for (int i = 0; i < diagonalMatrix.size(); i++)
	{
		ans[i][i] = diagonalMatrix[i];
	}

	return ans;
}

Matrix::Matrix()
{
}

Matrix::Matrix(int row, int col)
{
	_matrixContainer.resize(row, std::vector<double>(col, 0));
}

Matrix::Matrix(const MatrixContainerType &matrixContainer)
{
	_matrixContainer = matrixContainer;
}

Matrix::Matrix(const Matrix &matrix)
{
	_matrixContainer = matrix.GetContianer();
}

Matrix Matrix::operator*(const Matrix &matrix) const
{
	return Matrixf::MulMatrix(*this, matrix);
}

Matrix Matrix::operator+(const Matrix &matrix) const
{
	return Matrixf::AddMatrix(*this, matrix);
}

Matrix Matrix::operator-(const Matrix &matrix) const
{
	return Matrixf::SubMatrix(*this, matrix);
}

std::vector<double> &Matrix::operator[](int i)
{
	return _matrixContainer[i];
}

const std::vector<double> &Matrix::operator[](int i) const
{
	return _matrixContainer[i];
}

const Matrix::MatrixContainerType &Matrix::GetContianer() const
{
	return _matrixContainer;
}

int Matrix::Row() const
{
	return _matrixContainer.size();
}

int Matrix::Col() const
{
	return _matrixContainer[0].size();
}

Matrix &Matrix::Inverse()
{
	*this = Matrixf::InverseMatrix(this->_matrixContainer);
	return *this;
}

Matrix &Matrix::Transform()
{
	*this = Matrixf::TransposeMatrix(this->_matrixContainer);
	return *this;
}

Matrix::Matrix(const Matrixf::DiagonalMatrix &diagonalMatrix)
{
	*this = Matrixf::CreateDigonal(diagonalMatrix);
}

void Matrix::Show()
{
	for(int i = 0;i<Row();i++)
	{
		for(int j =0;j<Col();j++)
		{
			printf("%7.4f\t\t",(*this)[i][j]);
		}
		printf("\n");
	}
}

void Matrix::Normalize()
{
	for(int i = 0;i<Row();i++)
	{
		for(int j =0;j<Col();j++)
		{
			(*this)[i][j] = Mathf::Normalize((*this)[i][j]);
		}
	}
}